Equilibrium fluid-solid coexistence of hard spheres 
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We present a tethered Monte Carlo simulation of the crystallization of hard spheres. Our method 
boosts the traditional umbrella sampling to the point of making practical the study of constrained 
Gibbs' free energies depending on several crystalline order-parameters. We obtain high-accuracy 
estimates of the fluid-crystal coexistence pressure for up to 2916 particles (enough to accommodate 
fluid-solid interfaces). We are able to extrapolate to infinite volume the coexistence pressure (pco = 
11.5727(10)fcBr/(j^) and the interfacial free energy (7{ioo} = 0.636(ll)feBr/o-^). 



PACS numbers: 05.10.Ln,64.60.-i, 64.60.My, 64.70.D-. 

Crystallization is a vast field of research, where ex- 
periments and theory cross-fertilize. Hard-spheres (HS) 
provide a celebrated example: the numerical finding of a 
fluid-solid phase transition [1] motivated experiments on 
colloids [2] 13] • Finding an accurate procedure to locate 
the equilibrium phase boundaries for HS is a crucial step 
to address the self-assembly of complex molecules [4], 
as modeled by HS plus non-spherical interactions (e.g 
patchy [5| and Janus particles [S])- 

Up to now, numerical simulations of crystallization 
phase transitions have been well behind their fluid-fluid 
counterpart (e.g. vapor-liquid equilibria 7 ). Actually, 
HS are the preferred benchmark for numerical approaches 
to crystallization. Yet, the lack of exact solutions en- 
hances the importance of accurate numerical and/or ex- 
perimental studies. 

However, for preexisting numerical methods, a simula- 
tion whose starting configuration is a fluid never reaches 
the equilibrium crystal. Much as in experiments [3], the 
simulation gets stuck in a metastable crystal, or a de- 
fective crystal (or even a glass [S]). The proliferation of 
metastable states defeats optimized Monte Carlo (MC) 
methods that overcome free-energy barriers in simpler 
systems (SHU]. Besides, experimental and numerical de- 
terminations of the interfacial free energy are plainly in- 
consistent (maybe due to a small electrical charge in the 
colloidal particles [H]). 

Since feasible numerical methods [13j could not form 
the correct crystalline phase spontaneously, choosing the 
starting particle configuration became an issue (e.g. crys- 
talline or a carefully crafted mixture of solid and fluid 
phases). Methods can be classified as equilibrium or 
nonequilibrium. In the phase switch MC one tries 
to achieve fluid-crystal equilibrium (only up to iV = 500 
HS [15 ). An alternative is the separate computation of 
the fluid and solid free energies, supplemented with the 
conditions of equal pressure, temperature and chemical 
potential. For the fluid's free energy, one resorts to ther- 
modynamic integration, while choices are available for 
the crystal (Wigner-Seitz [16] , Einstein crystal [III [18] , 
Einstein molecule |19|). The nonequilibrium direct coex- 



istence method [201 121] handles larger systems [22] . 

As for the accuracy, in equilibrium computations the 
coexistence pressure Pco was obtained with precisions 
of ~ 0.1% (at finite N). Yet, the N values that can 
be simulated are rather small. An N ^ oo extrapola- 
tion is mandatory, which degrades the final accuracy to 
~ 1% [lU llSl [19] (results are summarized in Table |l|. 
The situation improves by an order of magnitude for the 
direct-coexistence method. With the exception of [15| . 
the different estimations of pco are compatible, although 
with widely differing accuracies. 

The computation of the interfacial free energy, 7, is 
more involved, since the issue of spatially heterogeneous 
mixtures of ffuid and solid can no longer be skipped (as 
done in equilibrium computations of Pco)- Indeed, recent 
estimations are either precise but mutually incompati- 
ble [23l , or of lesser accuracy [25] . 

Here, we introduce a tethered MC [2S1 ]^7\ approach 
to HS crystallization. The correct crystal appears in our 
simulation by constraining the value of two order param- 
eters. At variance with preexisting methods, the crystal 
found is independent from the starting particle configura- 
tion. Tethered MC provides a major simplification for the 
standard umbrella sampling method [2H1 [23] : chemical- 
potential differences among fiuid and crystal are very pre- 
cisely computed from a thermodynamic integration. In 
fact, our method resembles studies of liquid-vapor equi- 
libria [30l [31] . We go continuously from the fluid to the 
crystal by varying a reaction coordinate that labels the 
intermediate states. Rather than particle density, our 
reaction coordinate is a blend of bond-orientational crys- 
tal order parameters with different symmetries [32'-'34j . 
Very accurate determinations of the coexistence pressure 
and the interfacial free energy follow. The number of HS 
ranges 108 < N = An^ < 4000, {n integer). Our largest 
systems do show the surface-driven geometric transitions 
characteristic of the asymptotic large N regime [551137| . 

We consider N hard spheres of diameter cr, at con- 
stant pressure p, in a cubic box with periodic boundary 
conditions. The equilibrium crystal is face-centered cu- 
bic (FCC) [3H]. With the shorthand R for the particle 
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positions, {ri}fLi, Gibbs free energy g{p,T) is given by 

„-NI3g(p,T) _ 



dVe-^P^ J dR H{R) , (1) 



(A: de Broglie thermal wavelength, /3 = l/{kBT) and 
H{R) = if any pair of spheres overlaps, or 1 otherwise). 

We loosely constraint the values of two global order 
parameters, Qq and C. The well-known Qg detects the 
spatially coherent alignment of nearest-neighbors bonds 
in a lattice [321 ESI ■ It is the / = 6 instance of 
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{Y™{fij): spherical harmonics; fij: unitary vector point- 
ing from particle i to particle j; Nb{i): number of neigh- 
bors of particle Qe is positive in a crystal, while 
it is negligible {Qe ^ 1/VN) in a fluid. Yet, Qq's ro- 
tational invariance is a nuisance: enforcing a large Qe 
causes a crystal grain to grow in the fluid, but its orien- 
tation in the simulation box is arbitrary. In fact, when 
the grain finally hits itself through the periodic box's 
boundaries, long-lived metastable helicoidal crystals ap- 
pear. The cure is an order parameter with only cubic 
symmetry (Hj: 
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where c^{f) = x^y^{l- z^)+x^z^{l-y^)+y^z^{l-x^). 
C = 1 in an ideal, well aligned FCC, while C « for a 
fluid. Constraining a large C value suffices to obtain a 
nice crystal, irrespectively of the starting configuration 
(either a gas or an FCC structure). Still, helps us 
label unambiguously the intermediate states between the 
fluid and the FCC: some helicoidal crystals and the fluid- 
solid mixtures differ on their Qe values (but not on C). 

To enforce the quasi-constraints C{R) ~ C, Qe^R) « 
Qe [ISl m] , first multiply the integrand in Eq. ^ by 
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Na 



dQedCe" 



-[{Qe-Qe(R)Y + {C-C(R)Y] 



(4) 



The tunable parameter a tightens the quasi-constraints 
(we choose a = 200 [21] )• Exchanging the integration 
order in ([T]) yields 



-NI3g(p,T) 



dQe dCe- 



-NnM{Q<i,c,p) 



(5) 



where the effective potential^ QniQ^n C,p) is given by 
p(3Na 



dRdVuj{R,V;Q6,C,p), (6) 
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uj{R,V;Qe,C,p) being the tethered weight [1^ 

iO = H{R) e^'^P^"'^ [{Qe~Qe{R)f + {C~C(R)y] 




FIG. 1. (Color online) Vector field Vi?jv as computed from 
Eq. (|8|, for a system of A*' = 256 hard-spheres, at the coexis- 
tence pressure for the fluid-FCC phase transition (we scaled 
VJ7jv with a factor 1/a). Both the fluid and the FCC crystal 
are local minima of the effective potential, where VJ7jv = 0. 
The BCC coordinates are from N = 250. 



Our method relies on Fluctuation-Dissipation formu- 
lae flEl , obtained by taking derivatives in Eq. ([6| . We 
compute the gradient of at fixed pressure from: 

Vn^iQe, C) = a((Q6 - Qe{R)) , {C ~ C{R))) . (8) 

Coordinates (QgjC**) of local minima of f2 are lo- 
cated through 'VHn = 0. Furthermore, differences 
n]\[{Q^,C'') — ilNiQeiC"^) at fixed p are computed as 
the line integral of VJ7jv along any convenient path join- 
ing (0^,(7°) with iQlC'') in the (Qe,^) plane. 

The chemical potential g{p, T) is obtained from a 
saddle-point expansion in Eq. (|5|. Up to corrections 
vanishing as j3g{p,T) is the absolute minimum of 

^n{p,Q6,C). Yet, close to phase coexistence, fi^ has 
two relevant minima (i.e. the fluid and the FCC crys- 
tal). Therefore, the coexistence pressure pm^ follows 
from nf^'''^ 



O^'" (i.e. equal chemical potential). 



Our Metropolis MC simulation follows standard meth- 
ods [7]. We recast w in Eq. Q as the Boltzmann factor 
for HS at fixed pressure with a fictive potential energy 
kBTNa [(Qe - QeiR)? + (C - C{R)f]/2. Since Qe(-R) 
and C{R) are built out of sums of local terms, the num- 
ber of operations needed to compute their changes after 
a single-particle displacement does not grow with N . 

Our framework is illustrated in Fig. [T] where we show 

Vi7Ar(Qe, C) dX p = Pm \ We identify two local minima 



where V Qn = [the fluid, close to (Qe, C) = (1/V A^, 0), 
and the FCC minimum where both parameters are pos- 
itive]. Note their distance to other local minima of J7jv, 
such as the body centered cubic (BCC). 

Our main goal is to compute Ai7(p) — J?^*^^ — 
choosing the straight segment in Fig. [l] as integration 
path. The path is parameterized by our reaction coor- 
dinate, S {S^Q: fluid, S^^l: FCC). Actually, due to 
the additivity of Qe and C, choosing this segment is a 
must if we are to compute the interfacial free energy [41j . 
Indeed, physical fluid-solid coexistence is a convex com- 
bination of the two pure phases |42j, which provides a 
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FIG. 2. (Color online) (Top) Vi7jv projected over the liquid- 
FCC line, Vsi7jv, vs. the line parameter S {S — 0: fluid, S = 
1: FCC), for all our system sizes at the simulation pressures. 
(Bottom) Specific volume v = V/N as a. function of line 
parameter S. At large A'^, v becomes a linear function, as 
expected for a convex combination of pure phases |42) . 



physical interpretation for S as the fraction of particles 
in the coexisting solid phase: in the large N limit, v, C 
and Qe vary linearly with S (see Fig. [2] — bottom). 

Our simulation set up is as follows. We start by locat- 
ing {Qe, C) for the FCC and hquid minima at p ~ pio'^ ■ 
The first guess is obtained from NpT simulations with 
crystalline/disordered starting configurations. We later 
refine by solving for =0 [27 . 

Next, we introduce a uniform S grid on the liquid- 
FCC line and perform independent simulations at fixed 
{Qq,C,p) (see Appendix [a| for simulation details). As a 
test for equilibration, achieved for all N but N = 4000, 
every run was performed twice (starting from an ideal 
gas or from an ideal FCC crystal) [15] . 

Now, at variance with umbrella sampling, A^?(p) fol- 
lows from the integral over < S < 1 oi'V sf^N, the pro- 
jection of VJ7jv along the straight-line, Fig. [2] — top. We 
use reweighting extrapolations [271 144j to obtain AJ7(p) 

as a function of pressure. Then, it is easy to locate Pm\ 
Fig. |3] Statistical errors are estimated as in [TTj . 

We obtain pco = 11-5727(10) in units of k^T/a"^, in the 
large- TV limit. This result is six times more accurate than 
the best nonequilibrium estimate, Pco = 11.576(6) |22| 
and improves by a factor of 90 over the equilibrium esti- 
mate, Pco = 11.49(9) [13]. We compute Pco through a fair 
fit (x^ = 2.61 for three deg rees of freedom) of the ^co ^ 
listed in Table|l]to a second order polynomial in 1/7V [15] . 

As for the interfacial free energy, 7, we need to consider 
inhomogeneous configurations [15]. In fact, due to the 
periodic boundary conditions, at intermediate S the sur- 
face energy is minimized by mixed configurations where 
a crystalline slab (or cylinder, or bubble) is surrounded 
by fluid, see the snapshots in Appendix [B] As in vapor- 
liquid equilibria [311 137] , transitions among different ge- 
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FIG. 3. (Color online) Effective-potential difference Af2(p) = 
/jFCC __ ^fluid^ ^ function of pressure. At pi'^\ AQn = 0. 
The large iV limit stems from A^2(p) = [v^^^ - v^''''^){p - 
Pco)/{kBT) + 0{{p — Pco)^)- The simulated pressures (see 
Table [l]| correspond to the larger, filled symbols. 



ometries arise when S is varied. These transitions result 
in the cusps and steps that appear for large N in 'Vsf^N, 
Fig. [2] — top, and can be detected as well through the fluc- 
tuations of the particle density [27] . Under these circum- 
stances, 7 may be computed using Binder's method [47] . 
The effective potential has a local maximum along the 
line that joins the FCC and the fluid (the solution of 
Vs^N = at S* « 0.5, Fig. [2]— top). The excess 
free energy is due to the two interfaces that the fluid 
presents with a crystalline slab parallel to the simulation 
box ({100} planes). Then the interfacial free energy at 

(N) . 
Pco IS 

7^^) = fceT N (a* - nFCc) /(2 {Nv)%') . (9) 
The 7^^) (listed in Table |l| are extrapolated as [48] 
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(10) 



A fit for 256 < TV < 2916 yields 7 = 0.636(11) in units 
of k-oT/a^ (x^ = 0.14 for two degrees of freedom). We 
remark that the difference among the fit and 7(^=4000) 
is one fifth of the error bar (Table |l|. Also, the extrapo- 
lation for 500 < < 2916 merely doubles the final error 
estimate. Our result is compatible with 7 = 0.64(2) [5S] 
and 7 = 0.619(3) [M], but not with 7 = 0.5820(19) [33]. 
We remark that the 7^^^ estimation is fairly sensitive to 
p 27], an effect not systematically considered in [35H25] . 

Note that Eq. (10 1 holds if 7^^-' is computed at pco^-*. 

A final warning is in order. Not much is known about 
the effect of the Vs^^at's cusps and steps. Fig. [2]— top. 



in the large-TV extrapolation 7 
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This non- 



smoothness is a consequence of the geometric transitions 
that arise in our larger systems. However, the anal- 
ogy with simpler models [H] (e.g. the D = 2 Potts 
model, where comparison with exact solutions is possi- 
ble), strongly suggests that these cusps and steps are 
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Pco extrapolation. 
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TABLE I. For each TV, we report the simulated pressure in units of ksT/a^, the specific volume of the coexisting phases, the 
{100} surface tension 7{ioo} (in ksT/a^ units) and the phase-coexistence pressure (which is compared with work by other 
authors using different methods: phase switch Monte Carlo, the non-equilibrium direct coexistence method, and the Einstein 
Molecule approach). We extrapolate pi^' to the large iV limit as pi^'' = + ^i/N + Ui/N^ for 256 <A''< 2916 (p^"*°°° 
compatible but not included in the fit because of dubious equilibration) 
(Af>256 for the FCC and > 500 for the fluid). 



= p^ + ai/l\ + a2ll\' tor :^5b < JV < 29ie (,PcL"' — is 
The specific volume was extrapolated linearly in 1/N 



In summary, we have introduced a tethered MC p51[?7] 
approach to HS crystallization. We go continuously from 
the fluid to the crystal by varying a reaction coordi- 
nate. Tethered MC provides a major simplification to 
umbrella sampling, which makes it possible to study 
multi-constrained free energies. At variance with previ- 
ous methods, our simulations equilibrate (i.e. we find re- 
sults independent of the starting particle configuration) , 
not only for the formation of the space-filling crystal, 
but even for the more difficult case of mixed states with 
fluid-crystal interfaces. Our estimation of the coexistence 
pressure is, by far, the most accurate to date. That of 
the interfacial free energy is compatible with most (but 
not all) recent determinations. Should one wish to reach 
larger N , the tethered strategy would easily accommo- 
date additional order parameters. The method can also 
be generalized to other simple liquids, or to investigate 
the glass transition. 



As shown in the Fig. 2, we take a segment of the 
straight line in the (Qei C) plane that joins the fluid and 
the FCC minima of the effective potential. This segment 
is divided evenly in a grid of Ns points {Ns = 42 for 
N < 1372, while Ns = 82 for N > 2048). AU Ns points 
are simulated at the same pressure p (see Table 1^. We 
selected these p values by means of short, preliminary 
simulations. 

At each of the Ns {Qe,C,p) points, we run two in- 
dependent simulations with different initial conditions, 
an ideal gas and a perfect FCC crystal. Each of the 
2 X iVs runs had a length of 10'^ MC steps (1 MC step 
is composed of N particle displacements followed by one 
volume-change attempt). In addition, we show in Ta- 
ble |T] some A^— dependent observables computed in the 
simulation, as well as the large-A^ extrapolation. In or- 
der to ease comparison, we also tabulate the pca'^ values 
obtained by other groups (using different approaches). 
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Appendix A: Extended simulation details 

We provide some additional details for the interested 
reader. In particular, we give information necessary to 
reproduce our analysis and/or our simulations. 



Appendix B: Geometrical transitions 

As discussed in the main text, in a system with pe- 
riodic boundary conditions, geometrical transitions arise 
when the line parameter S varies from the liquid to the 
solid. In fact, the system struggles to minimize the sur- 
face energy while respecting the global constraints for 
Qq and C. Depending on the fraction of crystal phase, 
which is flxed by S', the minimizing geometry can be 
either a bubble, a cylinder or a slab of liquid in a crys- 
tal matrix (or vice versa). An example of each type of 
conflguration is displayed in Fig. |4] As S* varies, the 
minimizing geometry changes at definite S values. This 
phenomenon is named geometric transition^ and has been 
previously studied in simpler models (for instance, first- 
order transitions in lattice magnetic systems, or fluid-gas 
phase-coexistence) . 
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FIG. 4. Snapshots of mixed configurations for A'^ — 2916 particles found as the line parameter 5* varies. We present projections 
in the three Cartesian directions. To improve visibihty, the radii are a fraction of the real ones, and the darkness is an increasing 
function of the distance to the projection plane. 
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